A layered microcircuit model of somatosensory cortex with three interneuron types and cell-type-specific short-term plasticity

Abstract Three major types of GABAergic interneurons, parvalbumin-, somatostatin-, and vasoactive intestinal peptide-expressing (PV, SOM, VIP) cells, play critical but distinct roles in the cortical microcircuitry. Their specific electrophysiology and connectivity shape their inhibitory functions. To study the network dynamics and signal processing specific to these cell types in the cerebral cortex, we developed a multi-layer model incorporating biologically realistic interneuron parameters from rodent somatosensory cortex. The model is fitted to in vivo data on cell-type-specific population firing rates. With a protocol of cell-type-specific stimulation, network responses when activating different neuron types are examined. The model reproduces the experimentally observed inhibitory effects of PV and SOM cells and disinhibitory effect of VIP cells on excitatory cells. We further create a version of the model incorporating cell-type-specific short-term synaptic plasticity (STP). While the ongoing activity with and without STP is similar, STP modulates the responses of Exc, SOM, and VIP cells to cell-type-specific stimulation, presumably by changing the dominant inhibitory pathways. With slight adjustments, the model also reproduces sensory responses of specific interneuron types recorded in vivo. Our model provides predictions on network dynamics involving cell-type-specific short-term plasticity and can serve to explore the computational roles of inhibitory interneurons in sensory functions.


Introduction
Cortical GABAergic interneurons are inhibitory neurons that modulate and limit the degree of neuronal excitability in the neocortex.They can be classified according to electrophysiological or morphological characteristics or with molecular markers (Ascoli et al., 2008).Although challenges still exist in establishing consistency among different classification methods, many studies in recent years have used molecular markers, which label groups with different genetic origins, to make significant progress in exploring interneuron circuits (Tremblay et al., 2016;Campagnola et al., 2022).The three most common types of cortical interneurons express parvalbumin (PV), somatostatin (SOM), and vasoactive intestinal peptide (VIP), respectively, and play critical but distinct roles in cortical microcircuitry (Karnani et al., 2014;Tremblay et al., 2016;Campagnola et al., 2022).These three types have their own specific electrophysiology, morphology, connectivity, and short-term synaptic plasticity (STP).How these properties are related to their distinct dynamics and inhibitory functions is an important topic for understanding cortical microcircuitry.
PV and SOM cells have been extensively studied and compared experimentally.PV cells are considered a major stabilizing force that produces fast and reliable inhibition, mediated by synaptic activities that are weakened in amplitude (depressed) during high-frequency stimulation (Beierlein et al., 2003;Silberberg and Markram, 2007;Hu et al., 2014;Karnani et al., 2014;Naka and Adesnik, 2016).In contrast, SOM cells act more slowly, and have synaptic effects that are initially weak but can increase in amplitude (facilitate) in response to sustained high-frequency stimulation (Beierlein et al., 2003;Kapfer et al., 2007;Silberberg and Markram, 2007;Karnani et al., 2014;Yavorska and Wehr, 2016).Due to their different synaptic dynamics, PV and SOM cells contribute to cortical sensory processing (Natan et al., 2015(Natan et al., , 2017;;Seay et al., 2020) and control oscillatory activity (Chen et al., 2017;Van Derveer et al., 2021) in different but complementary ways.On the other hand, VIP cell activity results in disinhibition of pyramidal cells by inhibiting SOM cells (Lee et al., 2013;Pi et al., 2013;Karnani et al., 2016a).Distal or neuromodulatory inputs have been found to alter the activity of VIP cells to mediate disinhibition of sensory signals during different behavioral states, such as wakefulness and movements (Lee et al., 2013;Pi et al., 2013;Fu et al., 2014;Karnani et al., 2016a).Inhibitory interneurons thus form an integral part of the cortical computational circuitry, with different types contributing to different aspects of inhibitory control and complementing each other.
In addition, interneurons show a high degree of diversity across cortical layers.PV, SOM, and VIP cells have different morphologies and connectivities across layers (Xu et al., 2013;Prönneke et al., 2015;Tremblay et al., 2016;Feldmeyer et al., 2018).Furthermore, SOM cells exhibit different target preferences (Xu et al., 2013) and whisking-related activity (Muñoz et al., 2017) across layers of barrel cortex, resulting in layer-specific inhibitory modulation of network activity (Xu et al., 2013;Muñoz et al., 2017).The importance of this cross-layer diversity for computations in the cortical column remains largely unexplored (Hahn et al., 2022).
However, experiments are limited by the number of neurons that can be recorded in each animal or sample.In recent years, there has been a rapid development of genetic methods that allow different types of interneurons to be labeled and studied in living brains or brain slices (Josh Huang and Zeng, 2013;Huang, 2014).But even with genetic labeling, it is still relatively difficult to identify interneurons of a specific type in vivo, especially when it is located in deep cortical layers.Furthermore, in vivo paired recordings with specific interneuron types are even more technically challenging, which hinders the study of interneuron-typespecific synaptic dynamics.Therefore, to systematically study the roles of interneuron types, computational studies that take the interneuron diversity into account are needed in addition to experimental approaches.Computational models can help to develop a mechanistic understanding of interneuron functions and lead to hypotheses for experimental studies.In particular, understanding interneuron function and dynamics across layers is essential for the study of the multi-layered cortical microcircuitry.
A number of modeling studies have considered the roles of inhibitory neuron types in local cortical circuit dynamics.Yang et al. (2016) and Hertäg and Sprekeler (2019) examined how mutual inhibition between SOM and VIP cells allows switching between two processing modes in which top-down inputs to pyramidal cells are either integrated or canceled.Litwin- Kumar et al. (2016) studied PV, SOM, and VIP cells in firing rate models and revealed their dissociation in firing rate changes in the paradoxical effect, where stimulated inhibitory neurons show a paradoxical decrease in firing rate (Tsodyks et al., 1997;Murphy and Miller, 2009;Ozeki et al., 2009).Mahrach et al. (2020) compared experimental data of PV cell stimulation with their models and found that, compared to a simpler two-population (excitatory and inhibitory) model, a model with the three interneuron types can reflect more details in the experimental data, showing the paradoxical effect whether the network is inhibition-stabilized or not.Bos et al. (2020) further investigated the mechanics of a network with PV and SOM cells and showed how to predict their mean-field behaviors.Lee et al. (2017) and Lee and Mihalas (2017) studied the role of the three interneuron types in single-and multi-column sensory signal processing involving surround suppression and contexual modulation.While these studies have provided theoretical explanations for the functions of PV, SOM and VIP cells, a model that accounts for their diversity across layers is still needed to further understand their roles in a multi-layered cortical column.A few multi-layer models with interneuron types have been developed based on experimental data (Markram et al., 2015;Billeh et al., 2020;Borges et al., 2022;Moreni et al., 2023).However, incorporating the layerand cell-type-specific electrophysiological properties and synaptic dynamics of an interneuron type into the model while maintaining model simplicity and low simulation cost still remains a challenge, that, if achieved, will allow for more convenient use, adaptation, systematic analysis and generalization of the model.
We therefore developed a new cortical microcircuit model adapted from Potjans and Diesmann (2014), that incorporates layer-specific electrophysiological and synaptic properties of excitatory neurons (Exc) and the three interneuron types, PV, SOM and VIP.We focus on a mouse barrel column, taking parameters from mouse somatosensory cortex when available, complementing these with rat somatosensory cortex data.
To study the effects of cell-type-specific synaptic dynamics, we also include STP parameters derived from experimental data.Like Potjans and Diesmann (2014), we use a common and computationally low-cost neuron model, the Leaky Integrate-and-Fire (LIF) point neuron model.With the LIF neuron model we can already essentially capture resting-state spiking statistics and aspects of experimentally observed sensory responses.We further created a protocol of cell-type-specific stimulation in L2/3 and L4, by applying Poisson drives to all cells of a given type, to show network responses when different neuron types (Exc, PV, SOM, VIP) are activated.The results are consistent with the experimentally observed inhibitory effects of PV and SOM cells and the disinhibitory effects of VIP cells.In addition, the VIP population responses following stimulation of excitatory and the other inhibitory neurons provide predictions for future experiments.Using this protocol, we also compare results with and without STP, showing that STP may modulate population responses by altering the dominant inhibitory pathways.Specifically, the presence of STP qualitatively modifies the VIP cell responses to PV cell stimulation and the SOM cell responses to Exc cell stimulation.We hypothesize that this is because STP reduces the dominance of the PV→SOM→VIP and Exc→PV→SOM pathways.In summary, we created a simple, biologically plausible and computationally efficient model for the analysis of the roles of interneuron types in the barrel cortex and in multi-layered cortical microcircuitry more generally.

Model Overview
The model in this study is adapted from the multi-layer cortical microcircuit model by Potjans and Diesmann (2014).We extend the model to include three major interneuron types: the PV, SOM, and VIP cells, and use only experimental data of mouse and rat somatosensory cortex.Figure 1 shows an overview of the populations and their synaptic connectivity.All neurons are modeled as Leaky-Integrate-and-Fire (LIF) neurons with exponentially decaying postsynaptic currents (PSCs)1 (Tsodyks et al., 2000).Figure 2 shows the dimensions of the model, which are in accordance with those of a mouse C2 barrel column described previously (Lefort et al., 2009;Petersen, 2019).

Populations and Neuron Parameters
The model includes four cortical layers with Exc, PV, SOM, and VIP cells (Figure 1, Table 1).The cell numbers of populations in each layer are determined as follows.The layer-specific excitatory and inhibitory cell numbers (N Exc , N Inh ) are specified according to the estimation for the mouse C2 barrel column by Lefort et al. (2009).The PV, SOM, and VIP cell numbers (N PV , N SOM , N VIP ) in each layer are determined by distributing N Inh : where f PV , f SOM , f VIP are relative quantities of the three interneuron types, obtained from Figure 2D and 4D in Lee et al. (2010) with WebPlotDigitizer2 (see Relative Quantities of Interneuron Types, Supplementary Material).Since VIP cells are mostly located in L2/3 (Prönneke et al., 2015;Tremblay et al., 2016) and experimental data pertaining to their connectivity in other layers is lacking, we include all VIP cells in a single population for simplicity.This results in 13 populations in total (Table 1).
The cell-type-specific membrane parameters (Table 2) are based on the in vitro data by Neske et al. (2015).Their L2/3 data are used for L2/3 and L4 of our model, and their L5 data for L5 and L6, where we use weighted mean parameters for two excitatory subtypes according to their relative cell numbers.The membrane time constants (τ m ) are further adjusted to approximate the in vivo awake state as follows: Watson et al. (2008) reported a decrease in membrane resistance (R m ) by 50.9% for excitatory and by 4.9% for inhibitory neurons when they change from the Down state to the Up state.The Down and Up states are described as periods during which neuronal populations are silent and periods of long-duration multineuronal depolarization, respectively (Watson et al., 2008).We apply these reductions in R m to the τ m of the neurons in our model (Table 2), to approximate a change from in vitro to the in vivo awake state.
The absolute refractory period (τ ref ) is 2.0 ms for every neuron, as in Potjans and Diesmann (2014).(Petersen, 2019).The layer thicknesses are according to Lefort et al. (2009).Experimentally observed connection probabilities (Pexp) are adjusted to derive model connection probabilities (P , see Figure 3) that correspond to these dimensions (see Derivation of Connection Probabilities, Supplementary Data).

Synaptic Parameters
Synaptic transmission events in all connections are modeled as currents with an instantaneous rise and a monoexponential decay.Cortical EPSPs consist of mainly the AMPA and the NMDA receptor-mediated components (Feldmeyer et al., 2002); however, the AMPA component is the only one considered here, since we aim for modeling the resting state, during which many NMDA receptors are presumably blocked due to their voltage dependence in vivo (Mayer et al., 1984).The synaptic weights (w) are determined as follows.
First, the amplitudes of postsynaptic potentials (PSPs) are determined (Table 3).EPSPs of intracortical connections are firstly set to 0.5 mV, which is consistent with the range of in vivo recordings known to us (Jouhanneau et al., 2015(Jouhanneau et al., , 2018;;Pala and Petersen, 2015;Schoonover et al., 2014), and IPSPs are set to 2.0 mV, which is four times as strong as the EPSPs, as in Potjans and Diesmann (2014).With this as a start, adjustments are then made for certain intracortical connections (Table 3) to fine-tune the model (see Parameter Optimization, Methods).EPSPs of the thalamic input follow the values of thalamocortical connections in vivo reported by Bruno and Sakmann (2006) (Table 4), with the exception of SOM cells, for which we use 50% of the others to better reflect the reported weaker thalamocortical connections in this group of neurons (Ji et al., 2016).The intracortical EPSPs and IPSPs and thalamic EPSPs are log-normally distributed in the model to be consistent with data from in vivo (Jouhanneau et al., 2015(Jouhanneau et al., , 2018;;Pala and Petersen, 2015;Schoonover et al., 2014) and in vitro (Song et al., 2005) experiments.For intracortical EPSPs and IPSPs, the standard deviations are set to the same magnitude as the means (e.g., 0.5±0.5 mV), which is also consistent with the in vivo data (Jouhanneau et al., 2015(Jouhanneau et al., , 2018;;Pala and Petersen, 2015;Schoonover et al., 2014), where the standard deviations are 62 to 172% of the magnitude of the means.EPSP amplitudes of the background inputs are fixed to a value of 0.5 mV.Subsequently, the amplitudes of exponentially decaying postsynaptic currents (PSCs) required to produce such PSPs in the different neuron types are calculated (Rotter and Diesmann, 1999;Maksimov et al., 2018): (2) where C m is the membrane capacitance, τ m is the membrane time constant, τ syn is the decay time constant of PSC, and a stands for τsyn τm .C m and τ m depend on the postsynaptic population (Table 2), and τ syn is 2 ms for excitatory (τ syn,Exc ) and 4 ms for inhibitory (τ syn,Inh ) connections (Table 4), which are chosen to approximate in vitro data from rat (Feldmeyer et al., 2002) and mouse (Ma et al., 2012), respectively.The resulting PSCs are shown in Figure S1A and are used as the synaptic weights in the model code.
The synaptic delays d Exc , d Inh , d th (Table 4) follow in vivo data of excitatory and inhibitory intracortical connections (Jouhanneau et al., 2018) and thalamocortical connections (Bruno and Sakmann, 2006).They are log-normally distributed to be consistent with the experimental data as well.The synaptic delay of the Poisson background input (d bg ) is fixed to the same value as the integration resolution (0.1 ms) because it is irrelevant for network activity.Table 4. Synaptic time constants and delays.The time constants 1 are chosen to approximate in vitro data of rats and mice (Feldmeyer et al., 2002;Ma et al., 2012).The delays 2 are according to in vivo data of rats and mice (Bruno and Sakmann, 2006;Jouhanneau et al., 2018).For dExc and d Inh , values of pyramidal and PV cells in Jouhanneau et al. (2018) are used, respectively.
2 Latencies from the presynaptic spike to the start of postsynaptic current.

Probabilities of Intracortical Connections
The intracortical connections are created with pairwise Bernoulli trials with projection-specific probabilities (P ).Here, a projection stands for all connections from one neuron population to another (e.g., L4 E→L4 PV) (Senk et al., 2022).P of different projections are determined based on experimentally observed connection probabilities (P exp ) in paired recordings of mouse barrel or somatosensory cortex in vitro (Lefort et al., 2009;Packer and Yuste, 2011;Xu et al., 2013;Pala and Petersen, 2015;Walker et al., 2016;Karnani et al., 2016b;Jouhanneau et al., 2018;Kapfer et al., 2007;Ma et al., 2012;Scala et al., 2019;Hu et al., 2011;Galarreta and Hestrin, 2002;Nigro et al., 2018;Hilscher et al., 2017).Since every experiment has its own spatial range of recording, each P exp has to be adjusted to derive a P that corresponds to the dimensions our model (Figure 2).The steps of the derivation of P from P exp are described in Derivation of Connection Probabilities, Supplementary Material.We use an exponential decay for the connection probability (Packer and Yuste, 2011;Perin et al., 2011), simplifying the distance dependence that is in reality shaped in more detail by factors including the geometry and density of the neurites.The pairwise Bernoulli rule creates only a single synapse between a pair of neurons, according to the given connection probability.Thus, it lumps together potentially multiple synapses (multapses; for terminology, see Senk et al. (2022)) forming each connection in the brain into a single synapse.By choosing larger synaptic weights to reflect potential multapses, the only dynamical difference that remains compared to modeling each synapse separately is that the corresponding transmission delays become identical.In return, the computational efficiency of the simulations is increased.Since paired recordings simultaneously activate all synapses of the source neuron onto the target neuron, the synaptic weights derived from these experiments naturally fit into this modeling scheme.
Where experimental data are not available, we use assumed or estimated values for connection probability (Figure 3).For the case of intra-layer projections, the following rules are applied.1. Exc→PV equals PV→Exc, considering the observation of reciprocal coupling (Koelbl et al., 2013;Hu et al., 2014).2. L5 Exc→SOM, SOM→Exc and SOM→SOM use the averages of L2/3 and L4, due to the lack of more precise constraints.3. L6 uses the same values as L5, also due to the lack of more precise constraints, for all projections except Exc→Exc, where experimental data are available (Lefort et al., 2009).For the case of inter-layer projections, an approach based on physiologically determined connection probabilities is not possible, for lack of the corresponding data.Instead, we use the estimations of connection probabilities by Markram et al. (2015) 3 , which are based on morphological data, as P exp to derive P .Following Figure 2 and Table 1 in Markram et al. (2015), we take the pyramidal cells, star pyramidal cells, and spiny stellate cells as Exc cells, the large basket cells and normal basket cells as PV cells, the Martinotti cells as SOM cells, and the double bouquet cells and bipolar cells as VIP cells.Accordingly, for each layer the connection probabilities of morphological subtypes are combined by weighted averaging, taking into account their cell numbers in the database, to obtain the connection probabilities for Exc, PV, SOM, and VIP cells.Overall, 29 of 43 intra-layer connection probabilities and 12 of 126 inter-layer connection probabilities are directly derived from experimental data.For the others, it is still necessary to use assumptions and estimations as described above.

Background Input
The background input (Exc bg ) for each cell is a homogeneous Poisson spike input with a fixed EPSP amplitude of 0.5 mV and a constant but cell-type-specific firing rate (r bg ).The r bg for each cell type is optimized to obtain plausible population firing rates (see Parameter Optimization, Methods).The background input is always present in all simulations in this study.

Thalamic Input
To model a thalamic input (Exc th ), we estimate the cell number of a barreloid in the ventral posteromedial (VPM) nucleus of mouse corresponding to the C2 whisker.We derive a number of 115, by dividing the total cortical neuron number of our model (=6448) by the "Ratio S1/VPM" of rat C2 (=56) from Table 1 in Meyer et al. (2013).However, with firings according to touch-evoked responses of VPM (see text below), this number of cells produces much smaller cortical responses (data not shown) in our model than in vivo (Yu et al., 2019).Therefore, we double it to create an Exc th of 230 cells for larger cortical responses.We consider this adjustment biologically plausible, as neurons in the barrel cortex can respond to several adjacent whiskers with similar response latencies, i.e., they have multi-whisker receptive fields (Bruno and Simons, 2002;Le Cam et al., 2011;Staiger and Petersen, 2021).
The Exc th is connected to Exc and PV cells in all layers, while SOM cells are targeted only in L4, and VIP cells (located only in L2/3 in our model) are not targeted, following the pattern of cell-type-specific thalamocortical projections reported by Ji et al. (2016).The layer-specific connection probabilities (Figure 3) are specified according to the in vivo data from Constantinople and Bruno (2013).However, they did not report the connection probability for L2/3.Therefore, we estimate it as follows.We calculate the L2/3-to-L4 ratio of the number of VPM synapses per cortical neuron according to the data from Oberlaender et al. (2012).Then, assuming the connection probability for the pairwise Bernoulli connectivity is proportional to the average number of synapses per neuron (i.e., assuming the same multiplicity of synapses per pair of cortical and thalamic cells for the two layers), the L4 connection probability in Constantinople and Bruno (2013) is multiplied with this ratio to obtain an estimation of the L2/3 connection probability (Figure 3, Exc th →L2/3 Exc and Exc th →L2/3 PV).The synaptic weights are specified according to the in vivo data from Bruno and Sakmann (2006) (Table 3).We use these probabilities and weights for Exc and PV cells.For SOM cells, we use 50% of these values (Figure 3, Table 3) to reflect their weaker innervation percentages and strengths reported by Ji et al. (2016).
The thalamic stimulation is generated according to the touch-evoked responses of VPM from Figure 3F in Yu et al. (2019).The data are obtained with WebPlotDigitizer and offset along the y-axis to zero at stimulus onset (Figure 4, gray dots), and then fitted with a log-normal function to produce a time course of firing rate with a resolution of 0.1 ms (Figure 4, black line).This firing rate time course is generated in the form of inhomogeneous Poisson spike trains in all cells of Exc th to study cortical responses to a transient thalamic input, in separate simulations from those for resting-state activity.In each simulation, the thalamic stimulation is applied every 1 s and is repeated 10 times.

Parameter optimization
The firing rates (r bg ) of background input and a part of the synaptic weights described in the preceding are optimized to obtain plausible population firing rates and responses to the thalamic input.The adjustment of the background inputs is done in a cell-type-specific but layer-independent manner to limit the number of fitted parameters and thereby increase the robustness of the model.First, different values of r bg are tested for Exc, PV, SOM, and VIP cells to find a combination that results in the smallest deviation of the population firing rates from the in vivo data (Yu et al., 2019) in L2/3 and L4.The deviation is quantified by the root-mean-square percentage error of the population firing rates, i.e., for each population, the difference between the simulated and the in vivo data is computed as a fraction of the in vivo value, and the root mean square of the resulting values is taken.L5 and L6 are excluded at this step because they always result in large deviations (data not shown), possibly due to insufficient experiment-based connection probability data and consequent lesser reproduction of the in vivo activity.For PV and VIP cells, the r bg is further reduced slightly, to allow a stronger SOM cell activity overall.The resulting r bg for Exc, PV, SOM, and VIP cells are used throughout this study.For the model with STP (see Short-Term Synaptic Plasticity, Methods), a lower r bg for SOM cells is used to maintain plausible SOM cell firing rates, since the STP enhances their activity (data not shown).
Second, synaptic weights are adjusted (Table 3) to further reduce the deviation in firing rates from Yu et al. (2019) and to produce plausible responses to the thalamic input as in Yu et al. (2019) as well.Specifically, for L2/3, all intracortical excitatory synaptic weights to Exc and VIP cells are decreased by 50% to reduce their responses to the thalamic input, and the Exc→PV weight is increased by 50% to shift the PV-SOM balance towards PV.For L4, the Exc→PV weight is increased by 50% to also shift the PV-SOM balance towards PV.For L5 and L6, all intracortical excitatory synaptic weights to PV and SOM cells are decreased by 50% to prevent them from being overactive; this also allows higher Exc cell firing rates.Finally, PV→Exc and PV→SOM weights in L6 are increased by 50% to prevent the Exc and SOM cells in this layer from being overactive.

Short-Term Synaptic Plasticity
To study the effects of cell-type-specific short-term synaptic plasticity, we created another version of the model based on the original version (as described in the preceding) but incorporating STP synapses, and compared the two model versions in all analyses.The STP synapses are implemented as in the Tsodyks model (Tsodyks et al., 2000).This includes a parameter determining the dynamics of synaptic release probability (U), a facilitation time constant (F), and a depression time constant (D).
Values of U, F, and D for different projections are determined based on experimental data of unitary PSPs that demonstrate STP (Kapfer et al., 2007;Ma et al., 2012;Lefort and Petersen, 2017).For each projection (e.g., L4 Exc→L2/3 Exc), we created a pair of neurons connected with an STP synapse, and set the membrane parameters of the postsynaptic neuron to those of the corresponding population (Table 2).Then, we let the presynaptic neuron fire at the same constant rate as in the corresponding experiment.With this pair of neurons, we ran repeated simulations scanning U, F, and D (U: 0.05-1.0 in steps of 0.05; F, D: 0-1000 ms in steps of 20 ms) to determine the best-fit parameters, i.e., the set of U, F, and D that yields the smallest simulation-versus-experiment root-mean-sqaure error (RMSE) of normalized PSP amplitudes (examples shown in Figure 5A).The resulting set of U, F, and D is used for this particular projection.Note that Lefort and Petersen (2017) and Ma et al. (2012) subtracted overlapping components of preceding PSPs before calculating PSP amplitudes.Such a subtraction is also performed for the corresponding simulation data to ensure an accurate fit (e.g., L4 Exc→L2/3 Exc and L5 Exc→L6 Exc shown in Figure 5A).
With STP incorporated, the synaptic weights change during simulation.Therefore, to compare the two versions of the model under similar conditions, we scale the synaptic weights of the model with STP so that they evolve to a state which approximates the default, i.e., weights in the original model with static synapses (see Figure S1).To compute the scaling factors, we followed a numerical approach: For each projection, we created a pair of neurons connected with an STP synapse with its fitted STP parameters, and set the membrane parameters of the postsynaptic neuron to those of the corresponding population (Table 2).Then, we let the presynaptic neuron fire at the same rate as the corresponding population in the model with static synapses (Figure 6A2).Here, we used fixed inter-spike intervals, as this enabled a closer approximation of the resting state of the original model than fitting with Poisson statistics or spike trains taken directly from the model with static synapses.With the pair of neurons, we ran repeated 5-s simulations while scaling the initial synaptic weight, until the last synaptic weight recorded from the simulation deviated less than 0.1 pA from the default.The weight scaling factor thus obtained is used for this particular projection (Figure 5C).Note that because this approach relies on the recorded population firing rates in the model with static synapses, the obtained scaling factors always change with the random seed used for simulation; Figure 5C only shows the result of one instance of simulation.
For the thalamic input with STP, because the input itself is transient, the above weight scaling approach does not apply.Instead, the weight parameter is simply scaled so that the effective initial weight is equal to the default.This is done by setting w U as the weight parameter, so that the effective initial synaptic weight is U • w U = w, where U is the parameter determining the dynamics and initial value (as implemented in NEST) of synaptic release probability.Therefore, the weight scaling factor for the thalamic input is 1 U .

Network Resting State
We compared the resting state of our model, simulated using only the background (homogeneous Poisson) input, with two in vivo criteria: 1.Data on layer-and cell-type-specific firing rates by Yu et al. (2019), derived for excitatory cells and the three interneuron types.2. Criteria on asynchronous irregular (AI) activity by Maksimov et al. (2018).AI activity corresponds to a state of low synchrony and irregular neuronal firings, which exists in networks of LIF neurons with sparse connectivity and balanced inputs (Brunel, 2000) as well as in in vitro and in vivo samples (Steriade et al., 2001;Cruikshank et al., 2007;Destexhe et al., 2007;Johnson and Buonomano, 2007;Okun and Lampl, 2008;Compte et al., 2009;Tan and Wehr, 2009;Chauvette et al., 2011;Chen et al., 2012;Dehghani et al., 2016).Based on in vivo data, Maksimov et al. (2018) proposed criteria on cortical AI states for modeling studies.The study computes the pairwise spike count correlation and the coefficient of variation of inter-spike intervals (CV ISI) for in vivo data from rat and macaque (Watson et al., 2016;Chu et al., 2014).While the values were largely similar for the two species, we here consider the rat data.For both correlation and CV ISI, averages and ranges across 13 recording sessions from awake rat frontal cortex were calculated.For CV ISI, the average across Up states in anesthetized rat motor cortex is further provided.These results along with others were proposed as a set of in vivo criteria for evaluating cortical circuit models.We evaluate our model according to these criteria.For each layer, 200 neurons are randomly chosen regardless of cell type, to calculate the average pairwise spike count correlation.The same number of neurons are chosen separately to calculate the CV ISI, with different firing rate criteria (see below).Data from t = 10 s to t = 15 s of the simulations are used, and average results of 20 simulations are calculated to compare with the criteria.As in Maksimov et al. (2018), for pairwise spike count correlation, bin width is 10 ms and neurons with no firings during the sampled period are excluded before selection.For CV ISI, neurons with a firing rate lower than 1 spikes/s during the sampled period are excluded before selection.

Cell-Type-Specific Stimulation
To study how the network responds to activation of different cell types, additional Poisson inputs with a fixed EPSP amplitude of 0.5 mV are applied to different cell types in L2/3 and L4.The duration and interval of this input are both 1 s, and it is repeated 20 times and for nine levels of firing rate (up to 2000 spikes/s, depending on cell type; see Figure 7).The resultant population firing rates in the same layers are calculated for the later half of the simulation period (t = 500 ms to t = 1000 ms after stimulus onset), avoiding initial transient effects.This protocol is performed for Exc, SOM, PV, and VIP populations and for L2/3 and L4 separately in dedicated simulations other than the ones described in the preceding.
Hardware and Software Configurations for Simulations NEST 3.3 (Spreizer et al., 2022) is used for model implementation and simulation.All simulations are executed on a computing cluster with 48-thread compute nodes running at 2.5 GHz.The simulation resolution is 0.1 ms in all cases.For each simulation, 24 threads on one compute node were used, with OpenMP4 for parallel computing.For a resting-state simulation of 15 s biological time, both model versions take approximately 0.3 core-hours (∼2 s build phase and <45 s simulation phase on 24 threads).
For all analyses, 20 simulations with different random seeds for the Poisson input and the network connectivity are performed, and the first ten seconds of simulation are always excluded.

Parameter Fits for Short-Term Synaptic Plasticity
To determine cell-type-specific short-term plasticity parameters, we fitted the Tsodyks model to experimental STP data, i.e., unitary PSPs that demonstrate STP, as illustrated in Fig. 5A. Figure 5B shows the fitted STP parameters.Projections not shown in the figure remain static.For Exc→Exc projections (Figure 5B, bottom), experimental STP data from different layers are taken from Lefort and Petersen (2017).For projections involving interneurons (Figure 5B, top), experimental STP data are limited to L2/3 and L4.Kapfer et al. (2007) reported data involving L2/3 pyramidal cells, fast-spiking (FS) cells (taken as PV cells in our model), and SOM cells (pyramidal→FS, pyramidal→SOM).Ma et al. (2012) reported data involving L4 excitatory regular-spiking (RS) cells, FS cells, and SOM cells (FS→RS, FS→FS, FS→SOM, SOM→RS, SOM→PV).Since there is no data from the other layers, for these projections we use the fitted parameters regardless of the layer in our model.Karnani et al. (2016b) also report experimental STP data involving different cortical interneuron types, including VIP cells.However, adding this dataset resulted in an overactive resting state of the model (results not shown); therefore, it is not incorporated in this study.In addition, Hu and Agmon (2016) report experimental STP data on cell-type-specific thalamocortical projections.The averaged data on FS and SOM cells are listed in their Figure 1, but those on excitatory cells are not provided.Therefore, we digitized the data on FS and SOM cells (cell-attached data in the case of SOM cells; see their Figure 1R) for the STP of Exc th →PV and Exc th →SOM projections in our model.For Exc th →E, we use the STP of L4 Exc→L2/3 Exc in our model as a substitution (Figure 5B), which we likewise consider as a feedforward projection in sensory processing.
Figure 5C shows the weight scaling factors for the model with STP compared to that with static synapses.Figure S1 shows the resulting resting-state weights compared to the static version.No weight scaling is performed for projections that remain static or have a probability of zero.
Background inputs are optimized to obtain plausible population firing rates (see Parameter Optimization, Methods).For the model with static synapses, the optimized r bg are 5000, 6600, 2500, 3400 spikes/s for Exc, PV, SOM, VIP, respectively.For the model with STP, only the r bg for SOM is changed to 2250 spikes/s to maintain plausible SOM cell firing rates.

Modeled Resting State
The resting states of the model versions with static synapses and with STP are shown in Figure 6.The firing rates across layers and cell types are comparable to the non-whisking state data in Yu et al. (2019) (Figure 6, A2 and B2).Most populations have a simulation-versus-experiment difference in mean firing rate less than 50% of the in vivo data.For the model with static synapses, only L5 PV (+155%), L6 PV (+80%), and SOM (−55%) have a difference larger than 50%.With STP, only L2/3 VIP (−65%), L5 PV (+188%), and L6 PV (+94%) have a difference larger than 50%.Note that the error bars of the simulation data in Figure 6 indicate standard deviations across 20 simulation instances, while those of the experimental data indicate standard errors across neurons.The error bars are thus only a measure of variation and are not intended to compare the simulation and experimental data.
To compare the asynchronous irregular (AI) activity of the model with the in vivo condition, we use the criteria of pairwise spike count correlation and the coefficient of variation of the inter-spike intervals (CV ISI) established by Maksimov et al. (2018) (see Asynchronous Irregular State, Methods), which are based on both awake and anesthetized states.The correlations in different layers of the model are mostly in the range of the awake state (Figure 6, A3 and B3), with the exception of L5 in the model with STP.On the other hand, the CV ISI is always between the awake and anesthetized states (Figure 6, A4 and B4).For each layer, 200 neurons are randomly selected from neurons with firing rates > 0 spikes/s and regardless of cell type, resulting in 19,900 pairs of neurons.(A4) Coefficients of variation of the inter-spike interval distributions.For each layer, 200 neurons are randomly selected from neurons with firing rates ≥ 1 spike/s and regardless of cell type.In (A3) and (A4), dashed lines show the in vivo criteria as follows.awake max, awake min: the maximum and minimum in the means of 13 recording sessions in frontal cortices of awake rats.anesth.Up state: the mean of recorded Up states in motor cortex neurons of anesthetized rats.These criteria are established by Maksimov et al. (2018).( A1) to (A4) and ( B1) to (B4) show the data for the models with static synapses and with STP, respectively.

Network Responses to Cell-Type-Specific Stimulation
Figure 7 shows how cell-type-specific stimulation in L2/3 and L4 changes network activity in the same layer.The interneurons exhibit inhibitory or disinhibitory effects on Exc cells.Figure 7A displays the results for the model with static synapses.In L2/3, PV and SOM activation evoke inhibition while VIP activation evokes disinhibition: The Exc firing rate changes by −72.7% when the PV firing rate changes by +50.6% (with r stim = 2000 spikes/s) and by −10.8% when the SOM firing rate changes by +73.1% (with r stim = 50 spikes/s); it changes by +7.0% when the VIP firing rate changes by +75.9% (with r stim = 200 spikes/s).In L4, PV activation evokes inhibition while SOM activation evokes disinhibition: The Exc firing rate changes by −76.8% when the PV firing rate changes by +36.7% (with r stim = 1000 spikes/s) and and by +5.7% when the SOM firing rate changes by +163.5% (with r stim = 200 spikes/s).
The contrast between inhibitory L2/3 SOM cells and disinhibitory L4 SOM cells here is consistent with the finding in the experiment by Xu et al. (2013).They suggested that a difference in SOM cell connectivity between L2/3 and L4 contributed to this contrast, since they found that in L2/3 the SOM cells inhibit the pyramidal cells more strongly than they inhibit PV cells, while in L4 they instead inhibit the PV cells more strongly.To test this hypothesis, we use the model with static synapses to examine the results of L4 SOM activation with lower L4 SOM→PV connection probability (Figure 8).We find that L4 SOM cells become inhibitory to E cells when the L4 SOM→PV connection probability is reduced to 50% of the original value, in contrast to their original disinhibitory effects.This result is consistent with the hypothesis proposed by Xu et al. (2013).
In the model with STP, the interneurons exhibit the same inhibitory or disinhibitory effects in L2/3 and L4 (Figure 7B) as their counterparts in the model with static synapses.In L2/3, the Exc firing rate changes by −68.9% when the PV firing rate changes by +126.7% (with r stim = 2000 spikes/s) and by −4.0% when the SOM firing rate changes by +20.2% (with r stim = 50 spikes/s); it changes by +7.4% when the VIP firing rate changes by +77.4% (with r stim = 200 spikes/s).In L4, the Exc firing rate changes by −73.2% when the PV firing rate changes by +61.4% (with r stim = 1000 spikes/s) and by +2.3% when the SOM firing rate changes by +135.8% (with r stim = 200 spikes/s).
Although mostly similar, there are still qualitative differences in the results between the two versions of model.L2/3 VIP cells in the model with static synapses are activated in response to weak PV stimulation (250 spikes/s) but become suppressed with stronger stimulation, while those in the model with STP are only suppressed.L4 SOM cells in the model with static synapses are suppressed in response to Exc cell stimulation, while those in the model with STP keep approximately the same firing rate.rnorm: resulting population firing rate calculated from t = 500 ms to t = 1000 ms after stimulus onset, normalized to the data at rstim = 0.The shaded areas show standard deviations across the randomized instances.In cases where no shaded area is visible, the standard deviation is so small that the shading overlaps with the line.Right: responses with half (18.15%) of the original L4 SOM→PV connection probability.Dashed black lines are added to distinguish increasing (left) and decreasing (right) Exc firing rates.Notations are the same as in Figure 7.

Network Responses to Thalamic Stimulation
Figure 9 shows the peristimulus time histograms (PSTH) of all populations of cortical neurons in response to the thalamic stimulation, and compares them with the in vivo data from Yu et al. (2019).The responses have a pattern similar to the in vivo data: the PV cells have the fastest and strongest responses in all layers, followed by weaker responses of Exc and then SOM cells.Note that the stimulus onset in our model is the start of thalamic neuronal firings, instead of a whisker stimulation as in the experiment.Thus, the signal relay from whisker to thalamus is not considered.Possibly due to this, the peaks of cortical responses in our model occur earlier than in the in vivo data.In addition, in our model the SOM cell responses are obvious only in L2/3 and L4, and all responding populations show shorter activation than the in vivo data (Figure 3F in Yu et al. (2019)).
With STP, SOM cell responses in L2/3 and L4 are stronger (Figure 9).In particular, even with lower resting-state firing rates, L4 SOM cells in the model with STP show stronger responses than in the model with static synapses (Figure 9, A3 and B3).
Figure 10 shows the normalized probability densities of mean spike latencies (time from the stimulus onset to the first spike) of the excitatory neurons in response to the thalamic stimulation, and compares them with the in vivo data in Constantinople and Bruno (2013).The peaks of our data show the same order across layers as the in vivo data: the L4 peak is the highest and appears almost at the same time as the L5 peak, and they are followed by the L2/3 peak.However, as in the PSTH results, the peaks in our data occur earlier than the peaks in the in vivo data.The simulation-versus-experiment differences in peak time are −6.05ms, −5.57ms, and −6.59 ms for L2/3, L4, and L5 in the model with static synapses, and −5.9 ms, −5.49ms, and −6.72 ms for L2/3, L4, and L5 in the model with STP, respectively.Part of this difference is presumably again associated with the fact that the signal relay from whisker to thalamus is not considered.1F in Constantinople and Bruno (2013).In (A) and (B), data are from the same simulations as in Figure 9, the bin width is 0.5 ms, and the waveforms are smoothed by the Savitzky-Golay filter with a window of 10 points and an order of 3.

Discussion
We developed a computational model of a multi-layer cortical microcircuit incorporating three major inhibitory interneuron types, the parvalbumin-, somatostatin-and vasoactive intestinal peptide-expressing cells (PV, SOM and VIP cells).The model is constrained by biologically plausible parameters from mouse and rat somatosensory (S1) cortex on these three interneuron types.It shows plausible resting-state activities and responses to modeled sensory stimulation.It also provides predictions about network responses when different neuron types are stimulated, and about the effects of cell-type-specific STP in the inhibitory control of the network.Therefore, this model can help to theoretically and systematically study the microcircuit functions and mechanisms involving interneuron types across multiple layers of the somatosensory cortex.In this section, we discuss the links to previous studies, the limitations of the model, and potential future work.

Model Parameters
Compared to the anesthetized and in vitro conditions, the awake state is generally considered to have lower membrane resistance (R m ) and hence a shorter membrane time constant (τ m ), because of frequent synaptic activity (Destexhe et al., 2003).Data indeed show that R m of pyramidal neurons in awake mice is lower than in vitro (Petersen, 2017).Since in vivo datasets containing τ m , R m , resting and threshold potentials of all interneuron types (PV, SOM, VIP) are lacking, we used a set of in vitro data (Neske et al., 2015) and made an adjustment for the awake state, using experimental data on changes in R m following transitions from Down to Up states (see Populations and Neuron Parameters, Methods).In vivo electrophysiological experiments with the three types of interneurons could provide the data necessary for testing and refining the corresponding assumptions.
Theoretically, EPSP and IPSP driving forces in neurons fluctuate with synaptic activity and associated membrane potential changes.However, in vitro and in vivo studies suggest that the effects of excitatory inputs sum close to linearly at the soma, presumably because of the isolation of individual inputs by dendritic branches and spines (Leger et al., 2005;Araya et al., 2006).Since our neuron model is a point neuron, we only model the effects of synaptic inputs on the somatic membrane potential and not local effects on the dendrites; and we therefore use current-based synapses to capture the linear summation.Nevertheless, this may have led to the excessive inhibition generated by PV cells in our model upon receiving transient thalamic input, because the electrical driving force was not accounted for (see Thalamic Stimulation, Discussion).
Limited by the availability of experimental data, we have to use estimations and assumptions for certain model parameters.Because data on interlaminar connection probabilities involving interneurons are lacking, we use the algorithmic estimates by Markram et al. (2015) to supplement this part of our model.Markram et al. (2015) distinguished the neuron types by morphology and provided a table for the correlation with molecular markers (Figure 2 and Table 1, Markram et al. (2015)).As described in Methods, we calculate average connection probabilities accordingly for the projections in our model (e.g., large basket cells and nest basket cells express PV; their connection probability data were averaged and used for the PV cells in our model), although this mapping may be not very precise.With this approach, the resulting connection probabilities of interlaminar inhibitory projections are mostly <5 % (Figure 3).Morphological studies have suggested that, although some interneurons have axons that are mostly confined to their layers of origin, others still have significant interlaminar projections (Thomson and Bannister, 2003;Kumar and Ohana, 2008;Packer et al., 2013).Further experimental data on functional connectivity can verify if the connection probabilities in our model fairly represent these interlaminar projections and improve the estimates where necessary.
In addition, fast-spiking (FS) cells are taken as PV cells to obtain data on connection probability or STP in some cases, where molecular marking is not done in the experiments (Kapfer et al., 2007;Hu et al., 2011;Ma et al., 2012).We also use assumptions for some intralaminar connection probabilities, mostly for L5 and L6, where experimental data are not available (see Probabilities of Intracortical Connections, Methods).It should further be noted that data on connection probability obtained from paired recordings in brain slices, such as those included in this study, may suffer from underestimation because of truncation of connections (Lefort et al., 2009;Stepanyants et al., 2009).In this regard, more experimental data combining advanced neuron classification and connectivity reconstruction methods (Kebschull et al., 2016;Gouwens et al., 2020) may improve the accuracy of the connectivity in our model.
The firing rates of the background inputs (r bg ) are optimized for obtaining plausible resting-state population firing rates.We found that the resulting r bg shows a pattern of PV > Exc > VIP > SOM, which is also found in some experimental data of long-range inputs to these cell types in S1 in recent studies (Martinetti et al., 2022;Naskar et al., 2021).This indicates that our optimization is biologically plausible at least in this respect.This level of background input is also similar to that in the model of Potjans and Diesmann (2014).The effective strength of the external drive onto excitatory cells is 5000 spikes/s × 0.5 mV/spike = 2500 mV/s, similar to the one in Potjans and Diesmann (2014), which is 2000 × 8 spikes/s × 0.15 mV/spike = 2400 mV/s.This similarity also holds for the average onto the inhibitory cells (2091 vs. 2220 mV/s).Compared to Potjans and Diesmann (2014), we use a stronger unitary weight, which is based on data from paired recording experiments, but the summed strength of the external input is approximately conserved.Note, however, that the spatial extent of our model is 0.06 mm 2 as opposed to the 1 mm 2 of the Potjans and Diesmann (2014) model, so that our circuit includes a smaller percentage of the sending neurons, and a larger percentage of the input to the neurons is provided by the external drive.That the external drive is nevertheless close to that in their model appears mostly related to the different neuron parameters, which are cell-type-specific in our model but were taken to be cell-type-independent in the model of Potjans and Diesmann (2014).
Incorporating further interneuron diversity should be considered as a next step.In addition to connectivity, SOM and PV cells also have very diverse morphologies across layers (Muñoz et al., 2017;Feldmeyer et al., 2018;Gouwens et al., 2019).Although the vast majority of VIP cells is located in L2/3, smaller numbers still exist in the deeper layers which show different dendritic and axonal projection patterns (Prönneke et al., 2015).It was also revealed recently that a possible subgroup of PV cells can mediate a thalamus-driven disinhibition in L4 (Hua et al., 2022).How this diversity of interneurons may contribute to the inhibitory control and computation in the cortical column can be investigated in future by extending and refining our model.

Short-Term Synaptic Plasticity
Through fitting of data on post-synaptic potentials to the model of Tsodyks et al. (2000), we systematically determined cell-type-specific parameters of short-term plasticity (STP) that may be useful for future modeling studies.As explained in Methods, the synaptic weights of the model with STP are scaled so that its resting-state synaptic weights are close to those set for the original model with static synapses.Consequently, the population firing rates of the two versions of the model are similar, since the other parameters are the same except for a slight difference in the background input for SOM cells.As shown in Results, both versions of the model have plausible population firing rates.Therefore, the model with STP can represent the in vivo state under synaptic depression and facilitation.This allows computational studies of STP effects in vivo, for instance studying the role of STP in the effects of transient stimulation (see also Thalamic Stimulation, Discussion).Since the original synaptic weights before depression and facilitation (e.g., measured from a silent network state) in vivo are difficult to measure or estimate, the incorporation of STP with realistic parameters in a network model as we have done provides a tool complementing what is possible experimentally.
Similar to the model with static synapses, the asynchronous irregular (AI) activity in the model with STP is close to that observed in vivo, as assessed using criteria compiled by Maksimov et al. (2018).The main difference compared to those criteria is that the pairwise correlation of neuronal firings in L5 is slightly higher (Figure 6B3).This may be related to the experimental observations of deep layer synchronized oscillations, especially in L5 (Silva et al., 1991;Le Bon-Jego and Yuste, 2007;Berger TK and Markram, 2010).Although synchronized cortical oscillations have sometimes been considered to be generated by intrinsic cell membrane mechanisms (Silva et al., 1991;Le Bon-Jego and Yuste, 2007;Hayut et al., 2011), they have also been linked to the specific connectivity and synaptic properties of PV and SOM cells (Le Bon-Jego and Yuste, 2007;Fanselow et al., 2008;Chen et al., 2017;Funk et al., 2017;Veit et al., 2017;Domhof and Tiesinga, 2021), which are included in our model.Therefore, in addition to the AI activity in general, our model can also be used to study the network mechanisms behind deep-layer synchronized oscillations, including the contributions of specific connectivity patterns and synaptic properties of interneurons.

Cell-Type-Specific Stimulation
Simulated cell-type-specific stimulation in L2/3 shows that PV and SOM cells are inhibitory and VIP cells are disinhibitory (Figure 7), consistent with expectations based on the experimental literature (Beierlein et al., 2003;Kapfer et al., 2007;Silberberg and Markram, 2007;Hu et al., 2014;Karnani et al., 2014;Naka and Adesnik, 2016;Yavorska and Wehr, 2016;Lee et al., 2013;Pfeffer et al., 2013;Pi et al., 2013;Karnani et al., 2016a).Several results regarding VIP cells may be of interest for further study: 1.With Exc cell stimulation, the VIP cells are suppressed.Hypothetically, this may be due to the activation of SOM cells and a consequent SOM→VIP inhibition, as this SOM→VIP projection is supported by experimental data (Karnani et al., 2016b) and implemented in our model as well.2. VIP cells in the model with static synapses are activated by weak PV cell stimulation, whereas they are suppressed at all stimulation strengths in the model with STP.We hypothesize that this initial suppression is due to direct PV→VIP inhibition being weaker than the disinhibition of VIP cells through a PV→SOM→VIP pathway, while the depressing PV→SOM projection in the model with STP overrides this dominance of the indirect pathway.These results and hypotheses regarding VIP cells should be examined in further experimental or theoretical studies.
With L4 cell-type-specific stimulation in our model, we observe the following: 1. L4 SOM cells show a disinhibitory effect on L4 Exc cells, in contrast to L2/3 SOM cells, which show an inhibitory effect (Figure 7).We believe this reflects the higher SOM→PV connection probability in L4 than in L2/3 in our model (36.30% vs. 11.81%, Figure 3).This is consistent with the experiment by Xu et al. (2013), where a contrast of inhibition vs. disinhibition of pyramidal cells by L2/3 vs. L4 SOM cells is also observed and considered to be a result of the target preference of L4 SOM cells for PV cells (see Comparisons with Relevant Models, Discussion).This shows that our model can reflect layer-specific roles of SOM cells.2. L4 Exc cell stimulation suppresses L4 SOM interneurons in the model with static synapses; but in the model with STP, L4 SOM cells keep approximately the same firing rate throughout different levels of L4 Exc cell stimulation.The deactivation in the model with static synapses may be due to the fact that the Exc→SOM projection is dominated by the disynaptic Exc→PV→SOM pathway, while the facilitating Exc→SOM projection and depressing PV→SOM projection in the model with STP overrides this dominance.This is another topic for further experimental or theoretical studies.
In both model versions, the sensitivity of Exc, PV, SOM, and VIP cell activity to the stimulation of their own populations is very different.The induced firing rates in L2/3 are ordered as SOM>VIP>Exc>PV for most stimulus strengths (Fig. 7).To the best of our knowledge, the available experimental data do not allow for direct comparison with this result, for lack of analogous cell-type-specific stimulation experiments.Therefore, the result provides another prediction, which can be examined by future experiments distinguishing the given interneuron types.A possible cause for the observed difference in sensitivity is the cell-type specificity of the membrane time constants in L2/3, which also follow the order SOM>VIP>Exc>PV (Table 2).This is because a larger membrane time constant increases the area under the PSPs onto the cell and hence the probability of the cell being brought to fire.With further analytical methods, our model can help to predict other neuronal or circuit-level factors behind this result.

Thalamic Stimulation
To demonstrate the capability of our model to simulate sensory responses, we compared the network responses to thalamic stimulation with in vivo data (Yu et al., 2019;Constantinople and Bruno, 2013).Our model shows similar patterns of layer-and cell-type-specific responses (Figure 9, 10).In addition, the model with STP shows more obvious SOM cell responses in L4 (Figure 9, A3 vs. B3), even though their resting state firing rates are lower than those in the model with static synapses (Figure 6, A2 vs. B2).This shows a possible contribution of STP to the responsiveness of SOM cells to sensory signals.
Nevertheless, a common anomaly in both the peristimulus time histogram (PSTH) and spike latency data is the lack of prolonged responses (after ∼20 ms).There may be a few causes of this discrepancy with respect to the experimental data.The first is our use of current-based synapses, which neglects changes in the driving force given by the distance between the membrane potential and the synaptic reversal potential.When receiving the transient thalamic input, PV cells in our model are activated most strongly and hence will tend to produce the largest inhibition to all cell types.In reality, PV cells in vivo may be less strongly activated by thalamic input, because a higher membrane potential of the PV cells when activated could reduce the strength of the compound thalamocortical EPSC they receive.To test this possibility, we need a synapse model that takes changes in driving force into account, be it classical conductance-based synapses or an intermediate form between current-and conductance-based synapses that considers the near-linear summation of excitatory inputs at the soma.However, to the best of our knowledge, data on cell-type-specific synaptic conductance changes in vivo during such sensory responses are lacking.The second potential cause of the observed discrepancies is that higher-order or non-specific thalamic nuclei such as the posterior medial nucleus (POm) may also contribute to longer responses (Zhang and Bruno, 2019).These inputs have not been considered here as we only incorporate the ventral posteromedial nucleus (VPM).Parameters of higher-order thalamic inputs are required to test this possibility.
The in vivo thalamocortical activation of L4 and L5 barrel cortex neurons to tactile stimulation usually precedes that of L2/3, but the comparison between L4 and L5 has yielded ambiguous results (Armstrong-James et al., 1992;Constantinople and Bruno, 2013), which may reflect comparable strengths of direct thalamic inputs to L4 and L5, or a diversity of excitatory neurons in L5 (Constantinople and Bruno, 2013).Our results on mean spike latencies are consistent with these observations, with the L4 and L5 response latencies largely overlapping with each other and preceding L2/3 (Figure 10).The simulation-versus-experiment differences in peak time are in the range 5-7 ms.This can be partly explained by the signal relay from whisker to thalamus, which takes approximately 2.5-3 ms ms (Yu et al., 2016(Yu et al., , 2019) ) and is not considered in our model.Our model can be used to analyze the effects of the parameters already incorporated as well as new parameters (e.g., excitatory neuron types) on layer-specific response latencies and the cortical responses to sensory stimulation more generally.

Comparisons with Relevant Models
In recent years, increasing attention has been devoted to incorporating the major interneuron types in modeling studies to understand cortical microcircuit dynamics and signal processing (Yang et al., 2016;Del Molino et al., 2017;Lee and Mihalas, 2017;Lee et al., 2017;Hertäg and Sprekeler, 2019;Litwin-Kumar et al., 2016;Mahrach et al., 2020;Sanzeni et al., 2020;Hertäg and Clopath, 2022;Borges et al., 2022;Guo and Kumar, 2023;Wagatsuma et al., 2023;Moreni et al., 2023).In particular, several multi-layer models of S1 and visual (V1) cortices incorporating multiple interneuron types have been developed (Markram et al., 2015;Billeh et al., 2020;Borges et al., 2022;Moreni et al., 2023).Morphological or physiological data from S1 (Markram et al., 2015;Billeh et al., 2020;Borges et al., 2022) or V1 (Billeh et al., 2020;Moreni et al., 2023) were used to derive neuron and connectivity parameters and to establish models with Leaky-Integrate-and-Fire (LIF) (Markram et al., 2015;Billeh et al., 2020;Moreni et al., 2023) or multi-compartment (Markram et al., 2015;Billeh et al., 2020;Borges et al., 2022) neurons.With these models, network synchrony (Markram et al., 2015), oscillatory activity (Moreni et al., 2023), and selective sensory responses (Billeh et al., 2020) were studied.Our model uses the LIF neuron model and exclusively relies on S1 data to enhance its self-consistency.It also incorporates cell-type-specific membrane parameters, connection probabilities, and STP based on experimental data.With such a computationally low-cost but biologically plausible model, the results already resemble experimental observations and provide insight into the functional roles of interneuron diversity and cell-type-specific synaptic STP in the cortical microcircuitry.Mahrach et al. (2020) studied the paradoxical effect, where stimulation reduces rather than increases the firing rate of inhibitory cells (Tsodyks et al., 1997;Murphy and Miller, 2009;Ozeki et al., 2009).They stimulated PV cells in vivo in mouse anterior lateral motor cortex and barrel cortex and compared the results with those of a computational model.Their model is able to reproduce the paradoxical effect found in experimental data and provides predictions on the underlying parameter values.Specifically, for their "Model 1", the paradoxical effect of PV cells depends on J SV (J EE J VS − J ES J VE ), where J XY stands for the interaction strength from population Y to X (X, Y ∈ {E, S, V}; E, S, V stand for Exc, SOM, and VIP cells, respectively).In their theory, if J EE is small enough to make J SV (J EE J VS − J ES J VE ) negative, the PV cells should show the paradoxical effect.Here, we make a comparison between their Model 1 and L2/3 of our model (with static synapses), which both have E, PV, SOM, and VIP cells.With our original parameters, the paradoxical effect is absent (Figure 7).As an effort to eliminate differences between their model and ours that could block the paradoxical effect, we tested our L2/3 network 1. with J EE being zero or very weak (down to 1/128 of the original), which predicts a paradoxical effect in their model, 2. without the extra projections (VIP→Exc, VIP→PV, VIP→VIP) and layers (L4 to L6) that are not present in their model, and 3. with very weak stimulation strength for PV cells (down to 1 spike/s with a PSP amplitude of 0.5 mV), but still found no paradoxical effect (data not shown).A possibility that we have not examined is the network size.Mahrach et al. (2020) and Sanzeni et al. (2020) both indicated that a network size smaller than the ones they studied (76,800 neurons and a cortical surface area of 500+ µm diameter, respectively) might fail to show paradoxical effects.Since our model only has 6,448 cortical neurons and a theoretical surface area of ∼276 µm diameter, future studies with an up-scaled model are needed to test this possibility.Bos et al. (2020) used a spiking neuron model to analyze the influence of PV and SOM cells and showed that the role of SOM cells depends on two particular pathways: When the SOM→PV→Exc pathway dominates, SOM cells are disinhibitory, whereas when SOM→Exc and PV→PV dominates, SOM cells are inhibitory.Experimental results indicate that SOM cells can indeed be inhibitory or disinhibitory depending on the circuitry, as they show inhibitory effects in L2/3 but disinhibitory effects in L4 (Xu et al., 2013).We tested the pathway dependence of these effects in our model, considering the L4 SOM cells.With our original connectivity, the L2/3 SOM cells are inhibitory and L4 SOM cells are disinhibitory (Figure 7), which is consistent with Xu et al. (2013).By changing the L4 SOM→L4 PV connection probability, we also found that the L4 SOM cells can be either inhibitory or disinhibitory (Figure 8), reproducing the findings of Bos et al. (2020).These results are consistent with the theories and experimental observations that layerand population-specific connectivity contributes to the layer-specific roles of SOM cells (Xu et al., 2013;Muñoz et al., 2017;Bos et al., 2020).It should be noted that the layer-specific roles may ultimately depend on the neuronal subgroups involved, especially Martinotti vs. non-Martinotti SOM cells (Xu et al., 2013;Emmenegger et al., 2018;Muñoz et al., 2017).Such SOM cell subgroups are not explicitly considered in our model yet, but may have been already reflected in the layer-specific connectivity data we used (Xu et al., 2013;Nigro et al., 2018;Scala et al., 2019).Incorporating SOM cell subgroups in our model will allow further studies on this topic.

Outlook
Our model can be used as a convenient template for computational studies of complex mechanisms such as the effects of neuromodulators on sensory signal processing.Being based on the LIF neuron model, mean-field analyses can provide mechanistic explanations of the population-level dynamics, and make predictions of network activities with different parameters or inputs without running simulations (Layer et al., 2022).This can help to facilitate the exploration of the effects of cell and connection properties of interneurons, and examine the hypotheses we have proposed on the main pathways underlying the suppression and enhancement of cell-type-specific activity (see Cell-Type-Specific Stimulation, Discussion).In particular, a mean-field analysis incorporating synaptic STP (Romani et al., 2006) may help reveal the mechanisms underlying the association between STP and the roles of different types of interneurons.A further future step can be to use data from public anatomical databases of barrel cortex, such as the BarrelCortexInSilico (Udvary et al., 2022) 5 , that can help improve the model with more details on connectivity and other parameters.These further developments will allow more refined and systematic explorations of the roles played by different types of interneurons in cortical circuits.
1 Author contributions statement H.J. and S.v.A. conceived the study.H.J., G.Q. and D.F. contributed to the interpretation of the experimental data.H.J. created and implemented the model, performed the simulations, produced the figures, and wrote the first manuscript draft.H.J. analyzed and refined the model under supervision of S.v.A.All authors contributed to the final manuscript.

Figure 1 .
Figure 1.Model overview.(A) Populations and connectivities of the model.Exc: excitatory neurons.PV, SOM, VIP: parvalbumin-, somatostatin-, and vasoactive intestinal peptide-expressing inhibitory interneurons.Exc bg : background input; excitatory Poissonian input to each neuron with a fixed weight and constant but cell-type specific firing rates.Exc th : thalamic input; excitatory Poissonian input with a transient time course of firing rate, generated by thalamic neurons connected to the network with cell-type-specific connection probability.For details of Exc bg and Exc th , see Background Input and Thalamic Input, Methods.Interneurons in each layer are grouped by red boxes to show their average external connections (lines to and from the box) and specific internal connections (lines within the box).Only projections with a connection probability ≥ 4% in the model are shown (thin lines: 4-8%; thick lines: ≥ 8%).Note that in some cases this diagram may only partly reflect the connectivity in reality, due to limitations in the availability of experimental data (see Model Parameters, Discussion).(B) L2/3 interneurons and their projections as an example of the specific connectivity of interneurons.As in (A), thin and thick lines show connection probability of 4-8% and ≥ 8% in the model, respectively.(C) Connection rules for excitatory (Exc) and inhibitory (Inh) populations regardless of layers and interneuron types, following the graphical notation introduced in Senk et al. (2022).The three interneuron types are considered together as Inh in this panel.Symbol definitions for (C) are as follows.Solid and dashed lines: deterministic and probabilistic connections.δ: one-to-one connectivity.p: pairwise Bernoulli connectivity.A: autapses allowed.M : multapses not allowed.w ∼ LN : log-normally distributed synaptic weights.w: fixed synaptic weight.d ∼ LN : log-normally distributed synaptic delays.

Figure 2 .
Figure2.Dimensions of the microcircuit model in µm.These dimensions are based on the data of a mouse C2 barrel column.The surface area is equivalent to 200×300 µm(Petersen, 2019).The layer thicknesses are according toLefort et al. (2009).Experimentally observed connection probabilities (Pexp) are adjusted to derive model connection probabilities (P , see Figure3) that correspond to these dimensions (see Derivation of Connection Probabilities, Supplementary Data).

Figure 3 .
Figure 3. Connection probabilities (P ) in %. * indicate data derived from paired recording experiments.* * indicate estimations involving morphological data, based onMarkram et al. (2015).Otherwise, data are based on assumptions.For the approaches to obtain these data, see Probabilities of Intracortical Connections, Methods, and Thalamic Input, Methods.Exc th : Thalamic input.Blanks show where probability is zero.

Figure 4 .
Figure 4. Evolution of thalamic (VPM) firing rate in response to whisker touch.Gray dots: Experimental data digitized from Figure 3F in Yu et al. (2019).Black line: The time course fitted with a log-normal function.This fitted time course is used for the thalamic input (Exc th ) for the model (see Thalamic Inputs, Methods).

Figure 5 .
Figure 5. Short-term plasticity parameters.(A) Examples of STP parameter fits for four projections.U, F, and D are STP parameters.U: parameter for release probability.F: facilitation time constant.D: depression time constant.Blue lines: Simulated PSPs with the best-fit STP parameters (the values shown above).Red dots: Experimental data on PSP amplitudes.All simulated and experimental data are normalized to the amplitudes of the first PSP.For each projection, the best-fit STP parameters are determined with the simulation-versus-experiment RMSE of the PSP amplitudes (see Short-Term Synaptic Plasticity, Methods).If the experimental data include a subtraction of preceding PSPs, such a subtraction is also done for the simulated PSPs (e.g., the L4 Exc→L2/3 Exc and L5 Exc→L6 Exc projections shown here; gray lines show the simulated PSPs before subtraction).(B) Best-fit STP parameters of all projections with STP.Top: Exc→Inh, Inh→Inh, and thalamocortical (Exc th ) projections.Bottom: Layer-specific Exc→Exc projections.These best-fit STP parameters are used throughout this study for the version of the model with STP.Projections that remain static are not shown in the figure.(C) Weight scaling factors of different projections for the model with STP.The values shown are used to scale the weights of different projections, in order to yield a resting state approximating the model with static synapses (see Short-Term Synaptic Plasticity, Methods).Blanks: Connection probability is zero or the synaptic weight is static.

Figure 6 .
Figure 6.Resting state of the model.(A1) Raster plot of neuronal firings.(A2) Population firing rates: simulation (bars) vs. awake non-whisking state data from mice(Yu et al., 2019) (black markers).The error bars for the simulation data stand for standard deviations across 20 simulation instances, while those for the experimental data stand for standard errors across neurons.Color code as in (A1).(A3) Pairwise spike count correlations.For each layer, 200 neurons are randomly selected from neurons with firing rates > 0 spikes/s and regardless of cell type, resulting in 19,900 pairs of neurons.(A4) Coefficients of variation of the inter-spike interval distributions.For each layer, 200 neurons are randomly selected from neurons with firing rates ≥ 1 spike/s and regardless of cell type.In (A3) and (A4), dashed lines show the in vivo criteria as follows.awake max, awake min: the maximum and minimum in the means of 13 recording sessions in frontal cortices of awake rats.anesth.Up state: the mean of recorded Up states in motor cortex neurons of anesthetized rats.These criteria are established byMaksimov et al. (2018).(A1) to (A4) and (B1) to (B4) show the data for the models with static synapses and with STP, respectively.

Figure 7 .
Figure 7. L2/3 and L4 network responses to cell-type-specific stimulation.Each graph shows the result of stimulating one particular cell type (indicated by dashed lines) in one layer.(A) Model with static synapses.(B) Model with STP.rstim: rate of stimulation applied to each cell.For each cell type, results of rstim = 0 and eight other levels of rstim are determined, and the data include 20 randomized simulation instances.The stimulation consists of Poisson spikes with a duration and an interval of one second, and 20 repeats are performed for each level of rstim.rnorm:resulting population firing rate calculated from t = 500 ms to t = 1000 ms after stimulus onset, normalized to the data at rstim = 0.The shaded areas show standard deviations across the randomized instances.In cases where no shaded area is visible, the standard deviation is so small that the shading overlaps with the line.

Figure 8 .
Figure 8.Comparison of L4 network responses to activation of L4 SOM cells for different L4 SOM→L4 PV connection probabilities.Left: responses with the original L4 SOM→PV connection probability (36.3%).Right: responses with half (18.15%) of the original L4 SOM→PV connection probability.Dashed black lines are added to distinguish increasing (left) and decreasing (right) Exc firing rates.Notations are the same as in Figure 7.

Figure 9 .
Figure 9. Network responses in terms of PSTH upon thalamic stimulation.(A1) Raster plot of neuronal firings.(A2) Cell-type-specific PSTHs (cross-layer averages).(A3) Layer-and cell-type-specific PSTHs.The PSTH bin width is 1 ms.Crosses mark the peaks of the in vivo data, digitized from Figure 3 in Yu et al. (2019).The stimulus is generated as described in Thalamic Input, Methods.Data shown are the average of 20 randomized simulation instances, each with 10 repeats of stimulation.(A1) to (A3) and (B1) to (B3) show the data for the models with static synapses and with STP, respectively.

Figure 10 .
Figure 10.Normalized probability densities of mean spike latencies of excitatory neurons in response to thalamic stimulation.(A) Model with static synapses.(B) Model with STP.(C) In vivo data digitized from Figure1FinConstantinople and Bruno (2013).In (A) and (B), data are from the same simulations as in Figure9, the bin width is 0.5 ms, and the waveforms are smoothed by the Savitzky-Golay filter with a window of 10 points and an order of 3.

Table 1 .
Lee et al. (2010)09) numbers are based on the estimations of excitatory and inhibitory neuron numbers byLefort et al. (2009)and the relative quantities of PV, SOM and VIP cells byLee et al. (2010).For simplicity, all VIP cells are moved to L2/3.The numbers in parentheses show the VIP cell numbers before combining them into a single population.

Table 3 .
Postsynaptic potentials.These PSP amplitudes are defined for the model with static synapses.The PSCs (FigureS1, left) required to produce such PSPs are calculated with Equation (2).In the model with STP, the PSPs and PSCs are dynamic, except for the background input.